// 这个文件中包含向量版本的GEMM实现

#include "common.h"

#include <assert.h>
#include <math.h>
#include <riscv_vector.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>

#define __STR(x) #x
#define _STR(x) __STR(x)
#define ___xy(x, y) x##y
#define __xy(x, y) ___xy(x, y)
#define ___xyz(x, y, z) __xy(x, y##z)
#define __xyz(x, y, z) ___xyz(x, y, z)
#define __MX 8
#define _MX __xy(m, __MX)

#define vfloat32m1_t __xyz(register vfloat32, m1, _t)
#define vle32_v_f32m1 __xy(__riscv_vle32_v_f32, m1)
#define vfmadd_vf_f32m1 __xy(__riscv_vfmadd_vf_f32, m1)
#define vfmul_vf_f32m1 __xy(__riscv_vfmul_vf_f32, m1)
#define vse32_v_f32m1 __xy(__riscv_vse32_v_f32, m1)
#define vfadd_vv_f32m1 __xy(__riscv_vfadd_vv_f32, m1)

#define vfloat32m8_t __xyz(register vfloat32, m8, _t)
#define vle32_v_f32m8 __xy(__riscv_vle32_v_f32, m8)
#define vfmadd_vf_f32m8 __xy(__riscv_vfmadd_vf_f32, m8)
#define vfmul_vf_f32m8 __xy(__riscv_vfmul_vf_f32, m8)
#define vse32_v_f32m8 __xy(__riscv_vse32_v_f32, m8)
#define vfadd_vv_f32m8 __xy(__riscv_vfadd_vv_f32, m8)


/**
 * @brief 在8个向量寄存器为一组的情况下读取向量寄存器可以处理的最大元素个数。
 * @return 元素个数
 */
static inline unsigned long get_max_vlen_e32m8(void) {
  unsigned long vlen;
  asm volatile("vsetvli	%0,zero,e32,m8,ta,ma" : "=r"(vlen)::);
  return vlen;
}

/**
 * @brief 向量寄存器不分组的情况下读取向量寄存器可以处理的最大元素个数。
 * @return 元素个数
 */
static inline unsigned long get_max_vlen_e32m1(void) {
  unsigned long vlen;
  asm volatile("vsetvli	%0,zero,e32,m1,ta,ma" : "=r"(vlen)::);
  return vlen;
}

static int total_l1_cache_size = 1024 * 32 / sizeof(float); // l1 cache size = 32k
static int total_l2_cache_size = 1024 * 256 / sizeof(float); // l2 cache size = 256k

/**
 * @brief 算法 2：RVV向量优化版本GEMM。
 * @param[in]     M     矩阵A的行数
 * @param[in]     N     矩阵B的列数
 * @param[in]     K     矩阵A的列数和矩阵B的行数
 * @param[in]     ALPHA 标量系数
 * @param[in]     A     矩阵A的指针
 * @param[in]     B     矩阵B的指针
 * @param[in,out] C     输出矩阵C的指针
 * @return void
 */
void gemm_nn_rvv(int M, int N, int K, float ALPHA, const float *A, const float *B,
                 float *C) {
  // 计算每个向量操作可以处理的元素数
  size_t vl = get_max_vlen_e32m1();
  int tail = 0;
  int i, j, k;

  vfloat32m1_t vb;
  vfloat32m1_t vc;
  for (i = 0; i < M; ++i) {
    for (k = 0; k < K; ++k) {
      register float a_part = ALPHA * A[i * K + k];
      for (j = 0; j + vl <= N; j += vl) {
        vb = vle32_v_f32m1(&B[k * N + j], vl);
        vc = vle32_v_f32m1(&C[i * N + j], vl);
        vc = vfmadd_vf_f32m1(vb, a_part, vc, vl);
        vse32_v_f32m1(&C[i * N + j], vc, vl);
      }
      if (tail = N - j) {
        // 加载 B 的一部分到向量
        vb = vle32_v_f32m1(&B[k * N + j], tail);
        // 加载 C 的当前部分到向量
        vc = vle32_v_f32m1(&C[i * N + j], tail);
        // 执行向量乘法和加法
        vc = vfmadd_vf_f32m1(vb, a_part, vc, tail);
        // 将结果存回 C
        vse32_v_f32m1(&C[i * N + j], vc, tail);
      }
    }
  }
}

/**
 * @brief 步进读导致性能变差的GEMM版本。
 * @param[in]     M     矩阵A的行数
 * @param[in]     N     矩阵B的列数
 * @param[in]     K     矩阵A的列数和矩阵B的行数
 * @param[in]     ALPHA 标量系数
 * @param[in]     A     矩阵A的指针
 * @param[in]     B     矩阵B的指针
 * @param[in,out] C     输出矩阵C的指针
 * @return void
 */
void gemm_nn_rvv_k(int M, int N, int K, float ALPHA, const float *A, const float *B,
                   float *C) {
  // 计算每个向量操作可以处理的元素数
  size_t vl = get_max_vlen_e32m1();
  int tail = 0;
  int i, j, k;

  vfloat32m1_t vb;
  vfloat32m1_t vc;
  for (i = 0; i < M; ++i) {
    for (j = 0; j + vl <= N; j += vl) {
      vc = vle32_v_f32m1(&C[i * N + j], vl);
      for (k = 0; k < K; ++k) {
        register float a_part = ALPHA * A[i * K + k];
        vb = vle32_v_f32m1(&B[k * N + j], vl);
        vc = vfmadd_vf_f32m1(vb, a_part, vc, vl);
      }
      vse32_v_f32m1(&C[i * N + j], vb, vl);
    }
  }
}

/**
 * @brief 使用顺序读跨步写的方式将矩阵B打包到New_B中，寄存器不分组版本。
 * @param[in]     B     矩阵B的指针
 * @param[out]    New_B 输出矩阵New_B的指针
 * @param[in]     K     矩阵B的行数
 * @param[in]     N     矩阵B的列数
 * @param[in]     vl    向量长度
 * @return void
 */
static void pack_B_m1(const float *B, float *New_B, int K, int N, int vl) {
  int j, k, tail;

  for (k = 0; k < K; ++k) {
    for (j = 0; j + vl <= N; j += vl) {
      vfloat32m1_t vb = vle32_v_f32m1(&B[k * N + j], vl);
      vse32_v_f32m1(&New_B[j * K + vl * k], vb, vl);
    }
    if ((tail = N - j)) {
      vfloat32m1_t vb = vle32_v_f32m1(&B[k * N + j], tail);
      vse32_v_f32m1(&New_B[j * K + tail * k], vb, tail);
    }
  }
}

/**
 * @brief 算法 3：打包B矩阵版本GEMM。
 * @param[in]     M     矩阵A的行数
 * @param[in]     N     矩阵B的列数
 * @param[in]     K     矩阵A的列数和矩阵B的行数
 * @param[in]     ALPHA 标量系数
 * @param[in]     A     矩阵A的指针
 * @param[in]     B     矩阵B的指针
 * @param[in,out] C     输出矩阵C的指针
 * @return void
 */
void gemm_nn_rvv_pack(int M, int N, int K, float ALPHA, const float *A,
                      const float *B, float *C) {
  // 计算每个向量操作可以处理的元素数
  size_t vl = get_max_vlen_e32m1();
  int tail = 0;
  int i, j, k;

  // 打包 B 矩阵
  float *New_B = (float *)aligned_alloc(64, N * K * sizeof(float));
  pack_B_m1(B, New_B, K, N, vl);

  vfloat32m1_t vc;
  vfloat32m1_t vb;
  for (i = 0; i < M; ++i) {
    for (j = 0; j + vl <= N; j += vl) {
      vc = vle32_v_f32m1(&C[i * N + j], vl);
      for (k = 0; k < K; ++k) {
        register float a_part = ALPHA * A[i * K + k];
        vb = vle32_v_f32m1(&New_B[j * K + vl * k], vl);
        vc = vfmadd_vf_f32m1(vb, a_part, vc, vl);
      }
      vse32_v_f32m1(&C[i * N + j], vc, vl);
    }
    if ((tail = N - j)) {
      vc = vle32_v_f32m1(&C[i * N + j], tail);
      for (k = 0; k < K; ++k) {
        register float a_part = ALPHA * A[i * K + k];
        vb = vle32_v_f32m1(&New_B[j * K + tail * k], tail);
        vc = vfmadd_vf_f32m1(vb, a_part, vc, tail);
      }
      vse32_v_f32m1(&C[i * N + j], vc, tail);
    }
  }

  free(New_B);
}

/**
 * @brief 使用顺序读跨步写的方式将矩阵B打包到New_B中，8个向量寄存器为一组版本。
 * @param[in]     B     矩阵B的指针
 * @param[out]    New_B 输出矩阵New_B的指针
 * @param[in]     K     矩阵B的行数
 * @param[in]     N     矩阵B的列数
 * @param[in]     vl    向量长度
 * @return void
 */
static void pack_B_m8(const float *B, float *New_B, int K, int N, size_t vl) {
  int j, k, tail;
  for (k = 0; k < K; ++k) {
    j = 0;
    for (; j + 2 * vl <= N; j += 2 * vl) {
      vfloat32m8_t vb_00 = vle32_v_f32m8(&B[k * N + j], vl);
      vfloat32m8_t vb_01 = vle32_v_f32m8(&B[k * N + j + vl], vl);

      vse32_v_f32m8(&New_B[j * K + (2 * vl) * k], vb_00, vl);
      vse32_v_f32m8(&New_B[j * K + (2 * vl) * k + vl], vb_01, vl);
    }
    for (; j + vl <= N; j += vl) {
      vfloat32m8_t vb = vle32_v_f32m8(&B[k * N + j], vl);
      vse32_v_f32m8(&New_B[j * K + vl * k], vb, vl);
    }
    if ((tail = N - j)) {
      vfloat32m8_t vb = vle32_v_f32m8(&B[k * N + j], tail);
      vse32_v_f32m8(&New_B[j * K + tail * k], vb, tail);
    }
  }
}

/**
 * @brief 算法 4：向量寄存器分组版本GEMM。
 * @param[in]     M     矩阵A的行数
 * @param[in]     N     矩阵B的列数
 * @param[in]     K     矩阵A的列数和矩阵B的行数
 * @param[in]     ALPHA 标量系数
 * @param[in]     A     矩阵A的指针
 * @param[in]     B     矩阵B的指针
 * @param[in,out] C     输出矩阵C的指针
 * @return void
 */
void gemm_nn_rvv_m8(int M, int N, int K, float ALPHA, const float *A,
                    const float *B, float *C) {
  // 计算每个向量操作可以处理的元素数
  size_t vl = get_max_vlen_e32m8();
  int tail = 0;
  int i = 0, j, k;

  // 打包 B 矩阵
  float *New_B = (float *)aligned_alloc(64, N * K * sizeof(float));
  pack_B_m8(B, New_B, K, N, vl);
  for (; i < M; ++i) {
    j = 0;
    for (; j + 2 * vl <= N; j += 2 * vl) {
      vfloat32m8_t vc_00 = vle32_v_f32m8(&C[i * N + j], vl);
      vfloat32m8_t vc_01 = vle32_v_f32m8(&C[i * N + j + vl], vl);
      k = 0;
      for (; k + 3 < K; k += 4) {
        register float a_part = ALPHA * A[i * K + k];
        vfloat32m8_t vb_00;
        vfloat32m8_t vb_01;

        vb_00 = vle32_v_f32m8(&New_B[j * K + (2 * vl) * k], vl);
        vb_01 = vle32_v_f32m8(&New_B[j * K + (2 * vl) * k + vl], vl);
        vc_00 = vfmadd_vf_f32m8(vb_00, a_part, vc_00, vl);
        vc_01 = vfmadd_vf_f32m8(vb_01, a_part, vc_01, vl);

        a_part = ALPHA * A[i * K + (k + 1)];
        vb_00 = vle32_v_f32m8(&New_B[j * K + (2 * vl) * (k + 1)], vl);
        vb_01 = vle32_v_f32m8(&New_B[j * K + (2 * vl) * (k + 1) + vl], vl);
        vc_00 = vfmadd_vf_f32m8(vb_00, a_part, vc_00, vl);
        vc_01 = vfmadd_vf_f32m8(vb_01, a_part, vc_01, vl);

        a_part = ALPHA * A[i * K + (k + 2)];
        vb_00 = vle32_v_f32m8(&New_B[j * K + (2 * vl) * (k + 2)], vl);
        vb_01 = vle32_v_f32m8(&New_B[j * K + (2 * vl) * (k + 2) + vl], vl);
        vc_00 = vfmadd_vf_f32m8(vb_00, a_part, vc_00, vl);
        vc_01 = vfmadd_vf_f32m8(vb_01, a_part, vc_01, vl);

        a_part = ALPHA * A[i * K + (k + 3)];
        vb_00 = vle32_v_f32m8(&New_B[j * K + (2 * vl) * (k + 3)], vl);
        vb_01 = vle32_v_f32m8(&New_B[j * K + (2 * vl) * (k + 3) + vl], vl);
        vc_00 = vfmadd_vf_f32m8(vb_00, a_part, vc_00, vl);
        vc_01 = vfmadd_vf_f32m8(vb_01, a_part, vc_01, vl);
      }
      for (; k + 1 < K; k += 2) {
        register float a_part = ALPHA * A[i * K + k];
        vfloat32m8_t vb_00;
        vfloat32m8_t vb_01;

        vb_00 = vle32_v_f32m8(&New_B[j * K + (2 * vl) * k], vl);
        vb_01 = vle32_v_f32m8(&New_B[j * K + (2 * vl) * k + vl], vl);
        vc_00 = vfmadd_vf_f32m8(vb_00, a_part, vc_00, vl);
        vc_01 = vfmadd_vf_f32m8(vb_01, a_part, vc_01, vl);

        a_part = ALPHA * A[i * K + (k + 1)];
        vb_00 = vle32_v_f32m8(&New_B[j * K + (2 * vl) * (k + 1)], vl);
        vb_01 = vle32_v_f32m8(&New_B[j * K + (2 * vl) * (k + 1) + vl], vl);
        vc_00 = vfmadd_vf_f32m8(vb_00, a_part, vc_00, vl);
        vc_01 = vfmadd_vf_f32m8(vb_01, a_part, vc_01, vl);
      }
      for (; k < K; ++k) {
        register float a_part = ALPHA * A[i * K + k];
        vfloat32m8_t vb_00 = vle32_v_f32m8(&New_B[j * K + (2 * vl) * k], vl);
        vfloat32m8_t vb_01 =
            vle32_v_f32m8(&New_B[j * K + (2 * vl) * k + vl], vl);
        vc_00 = vfmadd_vf_f32m8(vb_00, a_part, vc_00, vl);
        vc_01 = vfmadd_vf_f32m8(vb_01, a_part, vc_01, vl);
      }
      vse32_v_f32m8(&C[i * N + j], vc_00, vl);
      vse32_v_f32m8(&C[i * N + j + vl], vc_01, vl);
    }
    for (; j + vl <= N; j += vl) {
      vfloat32m8_t vc = vle32_v_f32m8(&C[i * N + j], vl);
      k = 0;
      for (; k + 3 < K; k += 4) {
        register float a_part;
        vfloat32m8_t vb;

        a_part = ALPHA * A[i * K + k];
        vb = vle32_v_f32m8(&New_B[j * K + vl * k], vl);
        vc = vfmadd_vf_f32m8(vb, a_part, vc, vl);

        a_part = ALPHA * A[i * K + (k + 1)];
        vb = vle32_v_f32m8(&New_B[j * K + vl * (k + 1)], vl);
        vc = vfmadd_vf_f32m8(vb, a_part, vc, vl);

        a_part = ALPHA * A[i * K + (k + 2)];
        vb = vle32_v_f32m8(&New_B[j * K + vl * (k + 2)], vl);
        vc = vfmadd_vf_f32m8(vb, a_part, vc, vl);

        a_part = ALPHA * A[i * K + (k + 3)];
        vb = vle32_v_f32m8(&New_B[j * K + vl * (k + 3)], vl);
        vc = vfmadd_vf_f32m8(vb, a_part, vc, vl);
      }
      for (; k + 1 < K; k += 2) {
        register float a_part;
        vfloat32m8_t vb;

        a_part = ALPHA * A[i * K + k];
        vb = vle32_v_f32m8(&New_B[j * K + vl * k], vl);
        vc = vfmadd_vf_f32m8(vb, a_part, vc, vl);

        a_part = ALPHA * A[i * K + (k + 1)];
        vb = vle32_v_f32m8(&New_B[j * K + vl * (k + 1)], vl);
        vc = vfmadd_vf_f32m8(vb, a_part, vc, vl);
      }
      for (; k < K; ++k) {
        register float a_part = ALPHA * A[i * K + k];
        vfloat32m8_t vb = vle32_v_f32m8(&New_B[j * K + vl * k], vl);
        vc = vfmadd_vf_f32m8(vb, a_part, vc, vl);
      }
      vse32_v_f32m8(&C[i * N + j], vc, vl);
    }
    if ((tail = N - j)) {
      vfloat32m8_t vc = vle32_v_f32m8(&C[i * N + j], tail);
      k = 0;
      for (; k + 3 < K; k += 4) {
        register float a_part;
        vfloat32m8_t vb;

        a_part = ALPHA * A[i * K + k];
        vb = vle32_v_f32m8(&New_B[j * K + tail * k], tail);
        vc = vfmadd_vf_f32m8(vb, a_part, vc, tail);

        a_part = ALPHA * A[i * K + (k + 1)];
        vb = vle32_v_f32m8(&New_B[j * K + tail * (k + 1)], tail);
        vc = vfmadd_vf_f32m8(vb, a_part, vc, tail);

        a_part = ALPHA * A[i * K + (k + 2)];
        vb = vle32_v_f32m8(&New_B[j * K + tail * (k + 2)], tail);
        vc = vfmadd_vf_f32m8(vb, a_part, vc, tail);

        a_part = ALPHA * A[i * K + (k + 3)];
        vb = vle32_v_f32m8(&New_B[j * K + tail * (k + 3)], tail);
        vc = vfmadd_vf_f32m8(vb, a_part, vc, tail);
      }
      for (; k + 1 < K; k += 2) {
        register float a_part;
        vfloat32m8_t vb;

        a_part = ALPHA * A[i * K + k];
        vb = vle32_v_f32m8(&New_B[j * K + tail * k], tail);
        vc = vfmadd_vf_f32m8(vb, a_part, vc, tail);

        a_part = ALPHA * A[i * K + (k + 1)];
        vb = vle32_v_f32m8(&New_B[j * K + tail * (k + 1)], tail);
        vc = vfmadd_vf_f32m8(vb, a_part, vc, tail);
      }
      for (; k < K; ++k) {
        register float a_part = ALPHA * A[i * K + k];
        vfloat32m8_t vb = vle32_v_f32m8(&New_B[j * K + tail * k], tail);
        vc = vfmadd_vf_f32m8(vb, a_part, vc, tail);
      }
      vse32_v_f32m8(&C[i * N + j], vc, tail);
    }
  }

  free(New_B);
}

/**
 * @brief 使用顺序读跨步写的方式将矩阵B打包到New_B中，寄存器不分组但进行了循环展开。
 * @param[in]     B     矩阵B的指针
 * @param[out]    New_B 输出矩阵New_B的指针
 * @param[in]     K     矩阵B的行数
 * @param[in]     N     矩阵B的列数
 * @param[in]     vl    向量长度
 * @return void
 */
static void pack_B_loop_unroll(const float *B, float *New_B, int K, int N,
                               size_t vl) {
  int j, k, tail;

  for (k = 0; k < K; ++k) {
    j = 0;
    for (; j + 3 * vl <= N; j += 3 * vl) {
      vfloat32m1_t vb_00 = vle32_v_f32m1(&B[k * N + j], vl);
      vfloat32m1_t vb_01 = vle32_v_f32m1(&B[k * N + j + vl], vl);
      vfloat32m1_t vb_02 = vle32_v_f32m1(&B[k * N + j + 2 * vl], vl);

      vse32_v_f32m1(&New_B[j * K + (3 * vl) * k], vb_00, vl);
      vse32_v_f32m1(&New_B[j * K + (3 * vl) * k + vl], vb_01, vl);
      vse32_v_f32m1(&New_B[j * K + (3 * vl) * k + 2 * vl], vb_02, vl);
    }
    for (; j + 2 * vl <= N; j += 2 * vl) {
      vfloat32m1_t vb_00 = vle32_v_f32m1(&B[k * N + j], vl);
      vfloat32m1_t vb_01 = vle32_v_f32m1(&B[k * N + j + vl], vl);

      vse32_v_f32m1(&New_B[j * K + (2 * vl) * k], vb_00, vl);
      vse32_v_f32m1(&New_B[j * K + (2 * vl) * k + vl], vb_01, vl);
    }
    for (; j + vl <= N; j += vl) {
      vfloat32m1_t vb = vle32_v_f32m1(&B[k * N + j], vl);
      vse32_v_f32m1(&New_B[j * K + vl * k], vb, vl);
    }
    if ((tail = N - j)) {
      vfloat32m1_t vb = vle32_v_f32m1(&B[k * N + j], tail);
      vse32_v_f32m1(&New_B[j * K + tail * k], vb, tail);
    }
  }
}

/**
 * @brief 循环展开但寄存器不分组版本。
 * @param[in]     M     矩阵A的行数
 * @param[in]     N     矩阵B的列数
 * @param[in]     K     矩阵A的列数和矩阵B的行数
 * @param[in]     ALPHA 标量系数
 * @param[in]     A     矩阵A的指针
 * @param[in]     B     矩阵B的指针
 * @param[in,out] C     输出矩阵C的指针
 * @return void
 */
void gemm_nn_rvv_loop_unroll(int M, int N, int K, float ALPHA, const float *A,
                             const float *B, float *C) {
  // 计算每个向量操作可以处理的元素数
  size_t vl = get_max_vlen_e32m1();
  int tail = 0;
  int i = 0, j, k;

  // 打包 B 矩阵
  float *New_B = (float *)aligned_alloc(64, N * K * sizeof(float));
  pack_B_loop_unroll(B, New_B, K, N, vl);
  for (i; i + 7 < M; i += 8) {
    j = 0;
    for (; j + 3 * vl <= N; j += 3 * vl) {
      vfloat32m1_t vc_00 = vle32_v_f32m1(&C[i * N + j], vl);
      vfloat32m1_t vc_01 = vle32_v_f32m1(&C[i * N + j + vl], vl);
      vfloat32m1_t vc_02 = vle32_v_f32m1(&C[i * N + j + 2 * vl], vl);
      vfloat32m1_t vc_10 = vle32_v_f32m1(&C[(i + 1) * N + j], vl);
      vfloat32m1_t vc_11 = vle32_v_f32m1(&C[(i + 1) * N + j + vl], vl);
      vfloat32m1_t vc_12 = vle32_v_f32m1(&C[(i + 1) * N + j + 2 * vl], vl);
      vfloat32m1_t vc_20 = vle32_v_f32m1(&C[(i + 2) * N + j], vl);
      vfloat32m1_t vc_21 = vle32_v_f32m1(&C[(i + 2) * N + j + vl], vl);
      vfloat32m1_t vc_22 = vle32_v_f32m1(&C[(i + 2) * N + j + 2 * vl], vl);
      vfloat32m1_t vc_30 = vle32_v_f32m1(&C[(i + 3) * N + j], vl);
      vfloat32m1_t vc_31 = vle32_v_f32m1(&C[(i + 3) * N + j + vl], vl);
      vfloat32m1_t vc_32 = vle32_v_f32m1(&C[(i + 3) * N + j + 2 * vl], vl);
      vfloat32m1_t vc_40 = vle32_v_f32m1(&C[(i + 4) * N + j], vl);
      vfloat32m1_t vc_41 = vle32_v_f32m1(&C[(i + 4) * N + j + vl], vl);
      vfloat32m1_t vc_42 = vle32_v_f32m1(&C[(i + 4) * N + j + 2 * vl], vl);
      vfloat32m1_t vc_50 = vle32_v_f32m1(&C[(i + 5) * N + j], vl);
      vfloat32m1_t vc_51 = vle32_v_f32m1(&C[(i + 5) * N + j + vl], vl);
      vfloat32m1_t vc_52 = vle32_v_f32m1(&C[(i + 5) * N + j + 2 * vl], vl);
      vfloat32m1_t vc_60 = vle32_v_f32m1(&C[(i + 6) * N + j], vl);
      vfloat32m1_t vc_61 = vle32_v_f32m1(&C[(i + 6) * N + j + vl], vl);
      vfloat32m1_t vc_62 = vle32_v_f32m1(&C[(i + 6) * N + j + 2 * vl], vl);
      vfloat32m1_t vc_70 = vle32_v_f32m1(&C[(i + 7) * N + j], vl);
      vfloat32m1_t vc_71 = vle32_v_f32m1(&C[(i + 7) * N + j + vl], vl);
      vfloat32m1_t vc_72 = vle32_v_f32m1(&C[(i + 7) * N + j + 2 * vl], vl);
      for (k = 0; k < K; ++k) {
        register float a_part_0;
        register float a_part_1;
        register float a_part_2;
        register float a_part_3;
        register float a_part_4;
        register float a_part_5;
        register float a_part_6;
        register float a_part_7;

        vfloat32m1_t vb_00;
        vfloat32m1_t vb_01;
        vfloat32m1_t vb_02;

        a_part_0 = A[i * K + k];
        a_part_1 = A[(i + 1) * K + k];
        a_part_2 = A[(i + 2) * K + k];
        a_part_3 = A[(i + 3) * K + k];
        a_part_4 = A[(i + 4) * K + k];
        a_part_5 = A[(i + 5) * K + k];
        a_part_6 = A[(i + 6) * K + k];
        a_part_7 = A[(i + 7) * K + k];

        vb_00 = vle32_v_f32m1(&New_B[j * K + (3 * vl) * k], vl);
        vb_01 = vle32_v_f32m1(&New_B[j * K + (3 * vl) * k + vl], vl);
        vb_02 = vle32_v_f32m1(&New_B[j * K + (3 * vl) * k + 2 * vl], vl);

        vc_00 = vfmadd_vf_f32m1(vb_00, a_part_0, vc_00, vl);
        vc_01 = vfmadd_vf_f32m1(vb_01, a_part_0, vc_01, vl);
        vc_02 = vfmadd_vf_f32m1(vb_02, a_part_0, vc_02, vl);
        vc_10 = vfmadd_vf_f32m1(vb_00, a_part_1, vc_10, vl);
        vc_11 = vfmadd_vf_f32m1(vb_01, a_part_1, vc_11, vl);
        vc_12 = vfmadd_vf_f32m1(vb_02, a_part_1, vc_12, vl);
        vc_20 = vfmadd_vf_f32m1(vb_00, a_part_2, vc_20, vl);
        vc_21 = vfmadd_vf_f32m1(vb_01, a_part_2, vc_21, vl);
        vc_22 = vfmadd_vf_f32m1(vb_02, a_part_2, vc_22, vl);
        vc_30 = vfmadd_vf_f32m1(vb_00, a_part_3, vc_30, vl);
        vc_31 = vfmadd_vf_f32m1(vb_01, a_part_3, vc_31, vl);
        vc_32 = vfmadd_vf_f32m1(vb_02, a_part_3, vc_32, vl);
        vc_40 = vfmadd_vf_f32m1(vb_00, a_part_4, vc_40, vl);
        vc_41 = vfmadd_vf_f32m1(vb_01, a_part_4, vc_41, vl);
        vc_42 = vfmadd_vf_f32m1(vb_02, a_part_4, vc_42, vl);
        vc_50 = vfmadd_vf_f32m1(vb_00, a_part_5, vc_50, vl);
        vc_51 = vfmadd_vf_f32m1(vb_01, a_part_5, vc_51, vl);
        vc_52 = vfmadd_vf_f32m1(vb_02, a_part_5, vc_52, vl);
        vc_60 = vfmadd_vf_f32m1(vb_00, a_part_6, vc_60, vl);
        vc_61 = vfmadd_vf_f32m1(vb_01, a_part_6, vc_61, vl);
        vc_62 = vfmadd_vf_f32m1(vb_02, a_part_6, vc_62, vl);
        vc_70 = vfmadd_vf_f32m1(vb_00, a_part_7, vc_70, vl);
        vc_71 = vfmadd_vf_f32m1(vb_01, a_part_7, vc_71, vl);
        vc_72 = vfmadd_vf_f32m1(vb_02, a_part_7, vc_72, vl);
      }

      vc_00 = vfmul_vf_f32m1(vc_00, ALPHA, vl);
      vc_01 = vfmul_vf_f32m1(vc_01, ALPHA, vl);
      vc_02 = vfmul_vf_f32m1(vc_02, ALPHA, vl);
      vc_10 = vfmul_vf_f32m1(vc_10, ALPHA, vl);
      vc_11 = vfmul_vf_f32m1(vc_11, ALPHA, vl);
      vc_12 = vfmul_vf_f32m1(vc_12, ALPHA, vl);
      vc_20 = vfmul_vf_f32m1(vc_20, ALPHA, vl);
      vc_21 = vfmul_vf_f32m1(vc_21, ALPHA, vl);
      vc_22 = vfmul_vf_f32m1(vc_22, ALPHA, vl);
      vc_30 = vfmul_vf_f32m1(vc_30, ALPHA, vl);
      vc_31 = vfmul_vf_f32m1(vc_31, ALPHA, vl);
      vc_32 = vfmul_vf_f32m1(vc_32, ALPHA, vl);
      vc_40 = vfmul_vf_f32m1(vc_40, ALPHA, vl);
      vc_41 = vfmul_vf_f32m1(vc_41, ALPHA, vl);
      vc_42 = vfmul_vf_f32m1(vc_42, ALPHA, vl);
      vc_50 = vfmul_vf_f32m1(vc_50, ALPHA, vl);
      vc_51 = vfmul_vf_f32m1(vc_51, ALPHA, vl);
      vc_52 = vfmul_vf_f32m1(vc_52, ALPHA, vl);
      vc_60 = vfmul_vf_f32m1(vc_60, ALPHA, vl);
      vc_61 = vfmul_vf_f32m1(vc_61, ALPHA, vl);
      vc_62 = vfmul_vf_f32m1(vc_62, ALPHA, vl);
      vc_70 = vfmul_vf_f32m1(vc_70, ALPHA, vl);
      vc_71 = vfmul_vf_f32m1(vc_71, ALPHA, vl);
      vc_72 = vfmul_vf_f32m1(vc_72, ALPHA, vl);
      vse32_v_f32m1(&C[i * N + j], vc_00, vl);
      vse32_v_f32m1(&C[i * N + j + vl], vc_01, vl);
      vse32_v_f32m1(&C[i * N + j + 2 * vl], vc_02, vl);
      vse32_v_f32m1(&C[(i + 1) * N + j], vc_10, vl);
      vse32_v_f32m1(&C[(i + 1) * N + j + vl], vc_11, vl);
      vse32_v_f32m1(&C[(i + 1) * N + j + 2 * vl], vc_12, vl);
      vse32_v_f32m1(&C[(i + 2) * N + j], vc_20, vl);
      vse32_v_f32m1(&C[(i + 2) * N + j + vl], vc_21, vl);
      vse32_v_f32m1(&C[(i + 2) * N + j + 2 * vl], vc_22, vl);
      vse32_v_f32m1(&C[(i + 3) * N + j], vc_30, vl);
      vse32_v_f32m1(&C[(i + 3) * N + j + vl], vc_31, vl);
      vse32_v_f32m1(&C[(i + 3) * N + j + 2 * vl], vc_32, vl);
      vse32_v_f32m1(&C[(i + 4) * N + j], vc_40, vl);
      vse32_v_f32m1(&C[(i + 4) * N + j + vl], vc_41, vl);
      vse32_v_f32m1(&C[(i + 4) * N + j + 2 * vl], vc_42, vl);
      vse32_v_f32m1(&C[(i + 5) * N + j], vc_50, vl);
      vse32_v_f32m1(&C[(i + 5) * N + j + vl], vc_51, vl);
      vse32_v_f32m1(&C[(i + 5) * N + j + 2 * vl], vc_52, vl);
      vse32_v_f32m1(&C[(i + 6) * N + j], vc_60, vl);
      vse32_v_f32m1(&C[(i + 6) * N + j + vl], vc_61, vl);
      vse32_v_f32m1(&C[(i + 6) * N + j + 2 * vl], vc_62, vl);
      vse32_v_f32m1(&C[(i + 7) * N + j], vc_70, vl);
      vse32_v_f32m1(&C[(i + 7) * N + j + vl], vc_71, vl);
      vse32_v_f32m1(&C[(i + 7) * N + j + 2 * vl], vc_72, vl);
    }
    for (; j + 2 * vl <= N; j += 2 * vl) {
      vfloat32m1_t vc_00 = vle32_v_f32m1(&C[i * N + j], vl);
      vfloat32m1_t vc_01 = vle32_v_f32m1(&C[i * N + j + vl], vl);
      vfloat32m1_t vc_10 = vle32_v_f32m1(&C[(i + 1) * N + j], vl);
      vfloat32m1_t vc_11 = vle32_v_f32m1(&C[(i + 1) * N + j + vl], vl);
      vfloat32m1_t vc_20 = vle32_v_f32m1(&C[(i + 2) * N + j], vl);
      vfloat32m1_t vc_21 = vle32_v_f32m1(&C[(i + 2) * N + j + vl], vl);
      vfloat32m1_t vc_30 = vle32_v_f32m1(&C[(i + 3) * N + j], vl);
      vfloat32m1_t vc_31 = vle32_v_f32m1(&C[(i + 3) * N + j + vl], vl);
      vfloat32m1_t vc_40 = vle32_v_f32m1(&C[(i + 4) * N + j], vl);
      vfloat32m1_t vc_41 = vle32_v_f32m1(&C[(i + 4) * N + j + vl], vl);
      vfloat32m1_t vc_50 = vle32_v_f32m1(&C[(i + 5) * N + j], vl);
      vfloat32m1_t vc_51 = vle32_v_f32m1(&C[(i + 5) * N + j + vl], vl);
      vfloat32m1_t vc_60 = vle32_v_f32m1(&C[(i + 6) * N + j], vl);
      vfloat32m1_t vc_61 = vle32_v_f32m1(&C[(i + 6) * N + j + vl], vl);
      vfloat32m1_t vc_70 = vle32_v_f32m1(&C[(i + 7) * N + j], vl);
      vfloat32m1_t vc_71 = vle32_v_f32m1(&C[(i + 7) * N + j + vl], vl);

      for (k = 0; k < K; ++k) {
        register float a_part_0;
        register float a_part_1;
        register float a_part_2;
        register float a_part_3;
        register float a_part_4;
        register float a_part_5;
        register float a_part_6;
        register float a_part_7;

        a_part_0 = A[i * K + k];
        a_part_1 = A[(i + 1) * K + k];
        a_part_2 = A[(i + 2) * K + k];
        a_part_3 = A[(i + 3) * K + k];
        a_part_4 = A[(i + 4) * K + k];
        a_part_5 = A[(i + 5) * K + k];
        a_part_6 = A[(i + 6) * K + k];
        a_part_7 = A[(i + 7) * K + k];
        vfloat32m1_t vb_00;
        vfloat32m1_t vb_01;

        vb_00 = vle32_v_f32m1(&New_B[j * K + (2 * vl) * k], vl);
        vb_01 = vle32_v_f32m1(&New_B[j * K + (2 * vl) * k + vl], vl);
        vc_00 = vfmadd_vf_f32m1(vb_00, a_part_0, vc_00, vl);
        vc_01 = vfmadd_vf_f32m1(vb_01, a_part_0, vc_01, vl);
        vc_10 = vfmadd_vf_f32m1(vb_00, a_part_1, vc_10, vl);
        vc_11 = vfmadd_vf_f32m1(vb_01, a_part_1, vc_11, vl);
        vc_20 = vfmadd_vf_f32m1(vb_00, a_part_2, vc_20, vl);
        vc_21 = vfmadd_vf_f32m1(vb_01, a_part_2, vc_21, vl);
        vc_30 = vfmadd_vf_f32m1(vb_00, a_part_3, vc_30, vl);
        vc_31 = vfmadd_vf_f32m1(vb_01, a_part_3, vc_31, vl);
        vc_40 = vfmadd_vf_f32m1(vb_00, a_part_4, vc_40, vl);
        vc_41 = vfmadd_vf_f32m1(vb_01, a_part_4, vc_41, vl);
        vc_50 = vfmadd_vf_f32m1(vb_00, a_part_5, vc_50, vl);
        vc_51 = vfmadd_vf_f32m1(vb_01, a_part_5, vc_51, vl);
        vc_60 = vfmadd_vf_f32m1(vb_00, a_part_6, vc_60, vl);
        vc_61 = vfmadd_vf_f32m1(vb_01, a_part_6, vc_61, vl);
        vc_70 = vfmadd_vf_f32m1(vb_00, a_part_7, vc_70, vl);
        vc_71 = vfmadd_vf_f32m1(vb_01, a_part_7, vc_71, vl);
      }
      vc_00 = vfmul_vf_f32m1(vc_00, ALPHA, vl);
      vc_01 = vfmul_vf_f32m1(vc_01, ALPHA, vl);
      vc_10 = vfmul_vf_f32m1(vc_10, ALPHA, vl);
      vc_11 = vfmul_vf_f32m1(vc_11, ALPHA, vl);
      vc_20 = vfmul_vf_f32m1(vc_20, ALPHA, vl);
      vc_21 = vfmul_vf_f32m1(vc_21, ALPHA, vl);
      vc_30 = vfmul_vf_f32m1(vc_30, ALPHA, vl);
      vc_31 = vfmul_vf_f32m1(vc_31, ALPHA, vl);
      vc_40 = vfmul_vf_f32m1(vc_40, ALPHA, vl);
      vc_41 = vfmul_vf_f32m1(vc_41, ALPHA, vl);
      vc_50 = vfmul_vf_f32m1(vc_50, ALPHA, vl);
      vc_51 = vfmul_vf_f32m1(vc_51, ALPHA, vl);
      vc_60 = vfmul_vf_f32m1(vc_60, ALPHA, vl);
      vc_61 = vfmul_vf_f32m1(vc_61, ALPHA, vl);
      vc_70 = vfmul_vf_f32m1(vc_70, ALPHA, vl);
      vc_71 = vfmul_vf_f32m1(vc_71, ALPHA, vl);
      vse32_v_f32m1(&C[i * N + j], vc_00, vl);
      vse32_v_f32m1(&C[i * N + j + vl], vc_01, vl);
      vse32_v_f32m1(&C[(i + 1) * N + j], vc_10, vl);
      vse32_v_f32m1(&C[(i + 1) * N + j + vl], vc_11, vl);
      vse32_v_f32m1(&C[(i + 2) * N + j], vc_20, vl);
      vse32_v_f32m1(&C[(i + 2) * N + j + vl], vc_21, vl);
      vse32_v_f32m1(&C[(i + 3) * N + j], vc_30, vl);
      vse32_v_f32m1(&C[(i + 3) * N + j + vl], vc_31, vl);
      vse32_v_f32m1(&C[(i + 4) * N + j], vc_40, vl);
      vse32_v_f32m1(&C[(i + 4) * N + j + vl], vc_41, vl);
      vse32_v_f32m1(&C[(i + 5) * N + j], vc_50, vl);
      vse32_v_f32m1(&C[(i + 5) * N + j + vl], vc_51, vl);
      vse32_v_f32m1(&C[(i + 6) * N + j], vc_60, vl);
      vse32_v_f32m1(&C[(i + 6) * N + j + vl], vc_61, vl);
      vse32_v_f32m1(&C[(i + 7) * N + j], vc_70, vl);
      vse32_v_f32m1(&C[(i + 7) * N + j + vl], vc_71, vl);
    }
    for (; j + vl <= N; j += vl) {
      vfloat32m1_t vc_00 = vle32_v_f32m1(&C[i * N + j], vl);
      vfloat32m1_t vc_10 = vle32_v_f32m1(&C[(i + 1) * N + j], vl);
      vfloat32m1_t vc_20 = vle32_v_f32m1(&C[(i + 2) * N + j], vl);
      vfloat32m1_t vc_30 = vle32_v_f32m1(&C[(i + 3) * N + j], vl);
      vfloat32m1_t vc_40 = vle32_v_f32m1(&C[(i + 4) * N + j], vl);
      vfloat32m1_t vc_50 = vle32_v_f32m1(&C[(i + 5) * N + j], vl);
      vfloat32m1_t vc_60 = vle32_v_f32m1(&C[(i + 6) * N + j], vl);
      vfloat32m1_t vc_70 = vle32_v_f32m1(&C[(i + 7) * N + j], vl);

      for (k = 0; k < K; ++k) {
        register float a_part_0 = A[i * K + k];
        register float a_part_1 = A[(i + 1) * K + k];
        register float a_part_2 = A[(i + 2) * K + k];
        register float a_part_3 = A[(i + 3) * K + k];
        register float a_part_4 = A[(i + 4) * K + k];
        register float a_part_5 = A[(i + 5) * K + k];
        register float a_part_6 = A[(i + 6) * K + k];
        register float a_part_7 = A[(i + 7) * K + k];

        vfloat32m1_t vb_00 = vle32_v_f32m1(&New_B[j * K + vl * k], vl);
        vc_00 = vfmadd_vf_f32m1(vb_00, a_part_0, vc_00, vl);
        vc_10 = vfmadd_vf_f32m1(vb_00, a_part_1, vc_10, vl);
        vc_20 = vfmadd_vf_f32m1(vb_00, a_part_2, vc_20, vl);
        vc_30 = vfmadd_vf_f32m1(vb_00, a_part_3, vc_30, vl);
        vc_40 = vfmadd_vf_f32m1(vb_00, a_part_4, vc_40, vl);
        vc_50 = vfmadd_vf_f32m1(vb_00, a_part_5, vc_50, vl);
        vc_60 = vfmadd_vf_f32m1(vb_00, a_part_6, vc_60, vl);
        vc_70 = vfmadd_vf_f32m1(vb_00, a_part_7, vc_70, vl);
      }

      vc_00 = vfmul_vf_f32m1(vc_00, ALPHA, vl);
      vc_10 = vfmul_vf_f32m1(vc_10, ALPHA, vl);
      vc_20 = vfmul_vf_f32m1(vc_20, ALPHA, vl);
      vc_30 = vfmul_vf_f32m1(vc_30, ALPHA, vl);
      vc_40 = vfmul_vf_f32m1(vc_40, ALPHA, vl);
      vc_50 = vfmul_vf_f32m1(vc_50, ALPHA, vl);
      vc_60 = vfmul_vf_f32m1(vc_60, ALPHA, vl);
      vc_70 = vfmul_vf_f32m1(vc_70, ALPHA, vl);
      vse32_v_f32m1(&C[i * N + j], vc_00, vl);
      vse32_v_f32m1(&C[(i + 1) * N + j], vc_10, vl);
      vse32_v_f32m1(&C[(i + 2) * N + j], vc_20, vl);
      vse32_v_f32m1(&C[(i + 3) * N + j], vc_30, vl);
      vse32_v_f32m1(&C[(i + 4) * N + j], vc_40, vl);
      vse32_v_f32m1(&C[(i + 5) * N + j], vc_50, vl);
      vse32_v_f32m1(&C[(i + 6) * N + j], vc_60, vl);
      vse32_v_f32m1(&C[(i + 7) * N + j], vc_70, vl);
    }
    if ((tail = N - j)) {
      vfloat32m1_t vc_00 = vle32_v_f32m1(&C[i * N + j], tail);
      vfloat32m1_t vc_10 = vle32_v_f32m1(&C[(i + 1) * N + j], tail);
      vfloat32m1_t vc_20 = vle32_v_f32m1(&C[(i + 2) * N + j], tail);
      vfloat32m1_t vc_30 = vle32_v_f32m1(&C[(i + 3) * N + j], tail);
      vfloat32m1_t vc_40 = vle32_v_f32m1(&C[(i + 4) * N + j], tail);
      vfloat32m1_t vc_50 = vle32_v_f32m1(&C[(i + 5) * N + j], tail);
      vfloat32m1_t vc_60 = vle32_v_f32m1(&C[(i + 6) * N + j], tail);
      vfloat32m1_t vc_70 = vle32_v_f32m1(&C[(i + 7) * N + j], tail);

      for (k = 0; k < K; ++k) {
        register float a_part_0 = A[i * K + k];
        register float a_part_1 = A[(i + 1) * K + k];
        register float a_part_2 = A[(i + 2) * K + k];
        register float a_part_3 = A[(i + 3) * K + k];
        register float a_part_4 = A[(i + 4) * K + k];
        register float a_part_5 = A[(i + 5) * K + k];
        register float a_part_6 = A[(i + 6) * K + k];
        register float a_part_7 = A[(i + 7) * K + k];
        vfloat32m1_t vb_00 = vle32_v_f32m1(&New_B[j * K + tail * k], tail);
        vc_00 = vfmadd_vf_f32m1(vb_00, a_part_0, vc_00, tail);
        vc_10 = vfmadd_vf_f32m1(vb_00, a_part_1, vc_10, tail);
        vc_20 = vfmadd_vf_f32m1(vb_00, a_part_2, vc_20, tail);
        vc_30 = vfmadd_vf_f32m1(vb_00, a_part_3, vc_30, tail);
        vc_40 = vfmadd_vf_f32m1(vb_00, a_part_4, vc_40, tail);
        vc_50 = vfmadd_vf_f32m1(vb_00, a_part_5, vc_50, tail);
        vc_60 = vfmadd_vf_f32m1(vb_00, a_part_6, vc_60, tail);
        vc_70 = vfmadd_vf_f32m1(vb_00, a_part_7, vc_70, tail);
      }
      vc_00 = vfmul_vf_f32m1(vc_00, ALPHA, vl);
      vc_10 = vfmul_vf_f32m1(vc_10, ALPHA, vl);
      vc_20 = vfmul_vf_f32m1(vc_20, ALPHA, vl);
      vc_30 = vfmul_vf_f32m1(vc_30, ALPHA, vl);
      vc_40 = vfmul_vf_f32m1(vc_40, ALPHA, vl);
      vc_50 = vfmul_vf_f32m1(vc_50, ALPHA, vl);
      vc_60 = vfmul_vf_f32m1(vc_60, ALPHA, vl);
      vc_70 = vfmul_vf_f32m1(vc_70, ALPHA, vl);
      vse32_v_f32m1(&C[i * N + j], vc_00, tail);
      vse32_v_f32m1(&C[(i + 1) * N + j], vc_10, tail);
      vse32_v_f32m1(&C[(i + 2) * N + j], vc_20, tail);
      vse32_v_f32m1(&C[(i + 3) * N + j], vc_30, tail);
      vse32_v_f32m1(&C[(i + 4) * N + j], vc_40, tail);
      vse32_v_f32m1(&C[(i + 5) * N + j], vc_50, tail);
      vse32_v_f32m1(&C[(i + 6) * N + j], vc_60, tail);
      vse32_v_f32m1(&C[(i + 7) * N + j], vc_70, tail);
    }
  }
  for (i; i + 3 < M; i += 4) {
    j = 0;
    for (; j + 3 * vl <= N; j += 3 * vl) {
      vfloat32m1_t vc_00 = vle32_v_f32m1(&C[i * N + j], vl);
      vfloat32m1_t vc_01 = vle32_v_f32m1(&C[i * N + j + vl], vl);
      vfloat32m1_t vc_02 = vle32_v_f32m1(&C[i * N + j + 2 * vl], vl);
      vfloat32m1_t vc_10 = vle32_v_f32m1(&C[(i + 1) * N + j], vl);
      vfloat32m1_t vc_11 = vle32_v_f32m1(&C[(i + 1) * N + j + vl], vl);
      vfloat32m1_t vc_12 = vle32_v_f32m1(&C[(i + 1) * N + j + 2 * vl], vl);
      vfloat32m1_t vc_20 = vle32_v_f32m1(&C[(i + 2) * N + j], vl);
      vfloat32m1_t vc_21 = vle32_v_f32m1(&C[(i + 2) * N + j + vl], vl);
      vfloat32m1_t vc_22 = vle32_v_f32m1(&C[(i + 2) * N + j + 2 * vl], vl);
      vfloat32m1_t vc_30 = vle32_v_f32m1(&C[(i + 3) * N + j], vl);
      vfloat32m1_t vc_31 = vle32_v_f32m1(&C[(i + 3) * N + j + vl], vl);
      vfloat32m1_t vc_32 = vle32_v_f32m1(&C[(i + 3) * N + j + 2 * vl], vl);

      for (k = 0; k < K; ++k) {
        register float a_part_0 = A[i * K + k];
        register float a_part_1 = A[(i + 1) * K + k];
        register float a_part_2 = A[(i + 2) * K + k];
        register float a_part_3 = A[(i + 3) * K + k];
        vfloat32m1_t vb_00 = vle32_v_f32m1(&New_B[j * K + (3 * vl) * k], vl);
        vfloat32m1_t vb_01 =
            vle32_v_f32m1(&New_B[j * K + (3 * vl) * k + vl], vl);
        vfloat32m1_t vb_02 =
            vle32_v_f32m1(&New_B[j * K + (3 * vl) * k + 2 * vl], vl);
        vc_00 = vfmadd_vf_f32m1(vb_00, a_part_0, vc_00, vl);
        vc_01 = vfmadd_vf_f32m1(vb_01, a_part_0, vc_01, vl);
        vc_02 = vfmadd_vf_f32m1(vb_02, a_part_0, vc_02, vl);
        vc_10 = vfmadd_vf_f32m1(vb_00, a_part_1, vc_10, vl);
        vc_11 = vfmadd_vf_f32m1(vb_01, a_part_1, vc_11, vl);
        vc_12 = vfmadd_vf_f32m1(vb_02, a_part_1, vc_12, vl);
        vc_20 = vfmadd_vf_f32m1(vb_00, a_part_2, vc_20, vl);
        vc_21 = vfmadd_vf_f32m1(vb_01, a_part_2, vc_21, vl);
        vc_22 = vfmadd_vf_f32m1(vb_02, a_part_2, vc_22, vl);
        vc_30 = vfmadd_vf_f32m1(vb_00, a_part_3, vc_30, vl);
        vc_31 = vfmadd_vf_f32m1(vb_01, a_part_3, vc_31, vl);
        vc_32 = vfmadd_vf_f32m1(vb_02, a_part_3, vc_32, vl);
      }

      vc_00 = vfmul_vf_f32m1(vc_00, ALPHA, vl);
      vc_01 = vfmul_vf_f32m1(vc_01, ALPHA, vl);
      vc_02 = vfmul_vf_f32m1(vc_02, ALPHA, vl);
      vc_10 = vfmul_vf_f32m1(vc_10, ALPHA, vl);
      vc_11 = vfmul_vf_f32m1(vc_11, ALPHA, vl);
      vc_12 = vfmul_vf_f32m1(vc_12, ALPHA, vl);
      vc_20 = vfmul_vf_f32m1(vc_20, ALPHA, vl);
      vc_21 = vfmul_vf_f32m1(vc_21, ALPHA, vl);
      vc_22 = vfmul_vf_f32m1(vc_22, ALPHA, vl);
      vc_30 = vfmul_vf_f32m1(vc_30, ALPHA, vl);
      vc_31 = vfmul_vf_f32m1(vc_31, ALPHA, vl);
      vc_32 = vfmul_vf_f32m1(vc_32, ALPHA, vl);
      vse32_v_f32m1(&C[i * N + j], vc_00, vl);
      vse32_v_f32m1(&C[i * N + j + vl], vc_01, vl);
      vse32_v_f32m1(&C[i * N + j + 2 * vl], vc_02, vl);
      vse32_v_f32m1(&C[(i + 1) * N + j], vc_10, vl);
      vse32_v_f32m1(&C[(i + 1) * N + j + vl], vc_11, vl);
      vse32_v_f32m1(&C[(i + 1) * N + j + 2 * vl], vc_12, vl);
      vse32_v_f32m1(&C[(i + 2) * N + j], vc_20, vl);
      vse32_v_f32m1(&C[(i + 2) * N + j + vl], vc_21, vl);
      vse32_v_f32m1(&C[(i + 2) * N + j + 2 * vl], vc_22, vl);
      vse32_v_f32m1(&C[(i + 3) * N + j], vc_30, vl);
      vse32_v_f32m1(&C[(i + 3) * N + j + vl], vc_31, vl);
      vse32_v_f32m1(&C[(i + 3) * N + j + 2 * vl], vc_32, vl);
    }
    for (; j + 2 * vl <= N; j += 2 * vl) {
      vfloat32m1_t vc_00 = vle32_v_f32m1(&C[i * N + j], vl);
      vfloat32m1_t vc_01 = vle32_v_f32m1(&C[i * N + j + vl], vl);
      vfloat32m1_t vc_10 = vle32_v_f32m1(&C[(i + 1) * N + j], vl);
      vfloat32m1_t vc_11 = vle32_v_f32m1(&C[(i + 1) * N + j + vl], vl);
      vfloat32m1_t vc_20 = vle32_v_f32m1(&C[(i + 2) * N + j], vl);
      vfloat32m1_t vc_21 = vle32_v_f32m1(&C[(i + 2) * N + j + vl], vl);
      vfloat32m1_t vc_30 = vle32_v_f32m1(&C[(i + 3) * N + j], vl);
      vfloat32m1_t vc_31 = vle32_v_f32m1(&C[(i + 3) * N + j + vl], vl);

      for (k = 0; k < K; ++k) {
        register float a_part_0 = A[i * K + k];
        register float a_part_1 = A[(i + 1) * K + k];
        register float a_part_2 = A[(i + 2) * K + k];
        register float a_part_3 = A[(i + 3) * K + k];

        vfloat32m1_t vb_00 = vle32_v_f32m1(&New_B[j * K + (2 * vl) * k], vl);
        vfloat32m1_t vb_01 =
            vle32_v_f32m1(&New_B[j * K + (2 * vl) * k + vl], vl);
        vc_00 = vfmadd_vf_f32m1(vb_00, a_part_0, vc_00, vl);
        vc_01 = vfmadd_vf_f32m1(vb_01, a_part_0, vc_01, vl);
        vc_10 = vfmadd_vf_f32m1(vb_00, a_part_1, vc_10, vl);
        vc_11 = vfmadd_vf_f32m1(vb_01, a_part_1, vc_11, vl);
        vc_20 = vfmadd_vf_f32m1(vb_00, a_part_2, vc_20, vl);
        vc_21 = vfmadd_vf_f32m1(vb_01, a_part_2, vc_21, vl);
        vc_30 = vfmadd_vf_f32m1(vb_00, a_part_3, vc_30, vl);
        vc_31 = vfmadd_vf_f32m1(vb_01, a_part_3, vc_31, vl);
      }
      vc_00 = vfmul_vf_f32m1(vc_00, ALPHA, vl);
      vc_01 = vfmul_vf_f32m1(vc_01, ALPHA, vl);
      vc_10 = vfmul_vf_f32m1(vc_10, ALPHA, vl);
      vc_11 = vfmul_vf_f32m1(vc_11, ALPHA, vl);
      vc_20 = vfmul_vf_f32m1(vc_20, ALPHA, vl);
      vc_21 = vfmul_vf_f32m1(vc_21, ALPHA, vl);
      vc_30 = vfmul_vf_f32m1(vc_30, ALPHA, vl);
      vc_31 = vfmul_vf_f32m1(vc_31, ALPHA, vl);
      vse32_v_f32m1(&C[i * N + j], vc_00, vl);
      vse32_v_f32m1(&C[i * N + j + vl], vc_01, vl);
      vse32_v_f32m1(&C[(i + 1) * N + j], vc_10, vl);
      vse32_v_f32m1(&C[(i + 1) * N + j + vl], vc_11, vl);
      vse32_v_f32m1(&C[(i + 2) * N + j], vc_20, vl);
      vse32_v_f32m1(&C[(i + 2) * N + j + vl], vc_21, vl);
      vse32_v_f32m1(&C[(i + 3) * N + j], vc_30, vl);
      vse32_v_f32m1(&C[(i + 3) * N + j + vl], vc_31, vl);
    }
    for (; j + vl <= N; j += vl) {
      vfloat32m1_t vc_00 = vle32_v_f32m1(&C[i * N + j], vl);
      vfloat32m1_t vc_10 = vle32_v_f32m1(&C[(i + 1) * N + j], vl);
      vfloat32m1_t vc_20 = vle32_v_f32m1(&C[(i + 2) * N + j], vl);
      vfloat32m1_t vc_30 = vle32_v_f32m1(&C[(i + 3) * N + j], vl);

      for (k = 0; k < K; ++k) {
        register float a_part_0 = A[i * K + k];
        register float a_part_1 = A[(i + 1) * K + k];
        register float a_part_2 = A[(i + 2) * K + k];
        register float a_part_3 = A[(i + 3) * K + k];
        vfloat32m1_t vb_00 = vle32_v_f32m1(&New_B[j * K + vl * k], vl);
        vc_00 = vfmadd_vf_f32m1(vb_00, a_part_0, vc_00, vl);
        vc_10 = vfmadd_vf_f32m1(vb_00, a_part_1, vc_10, vl);
        vc_20 = vfmadd_vf_f32m1(vb_00, a_part_2, vc_20, vl);
        vc_30 = vfmadd_vf_f32m1(vb_00, a_part_3, vc_30, vl);
      }
      vc_00 = vfmul_vf_f32m1(vc_00, ALPHA, vl);
      vc_10 = vfmul_vf_f32m1(vc_10, ALPHA, vl);
      vc_20 = vfmul_vf_f32m1(vc_20, ALPHA, vl);
      vc_30 = vfmul_vf_f32m1(vc_30, ALPHA, vl);
      vse32_v_f32m1(&C[i * N + j], vc_00, vl);
      vse32_v_f32m1(&C[(i + 1) * N + j], vc_10, vl);
      vse32_v_f32m1(&C[(i + 2) * N + j], vc_20, vl);
      vse32_v_f32m1(&C[(i + 3) * N + j], vc_30, vl);
    }
    if ((tail = N - j)) {
      vfloat32m1_t vc_00 = vle32_v_f32m1(&C[i * N + j], tail);
      vfloat32m1_t vc_10 = vle32_v_f32m1(&C[(i + 1) * N + j], tail);
      vfloat32m1_t vc_20 = vle32_v_f32m1(&C[(i + 2) * N + j], tail);
      vfloat32m1_t vc_30 = vle32_v_f32m1(&C[(i + 3) * N + j], tail);

      for (k = 0; k < K; ++k) {
        register float a_part_0 = A[i * K + k];
        register float a_part_1 = A[(i + 1) * K + k];
        register float a_part_2 = A[(i + 2) * K + k];
        register float a_part_3 = A[(i + 3) * K + k];
        vfloat32m1_t vb_00 = vle32_v_f32m1(&New_B[j * K + tail * k], tail);
        vc_00 = vfmadd_vf_f32m1(vb_00, a_part_0, vc_00, tail);
        vc_10 = vfmadd_vf_f32m1(vb_00, a_part_1, vc_10, tail);
        vc_20 = vfmadd_vf_f32m1(vb_00, a_part_2, vc_20, tail);
        vc_30 = vfmadd_vf_f32m1(vb_00, a_part_3, vc_30, tail);
      }
      vc_00 = vfmul_vf_f32m1(vc_00, ALPHA, vl);
      vc_10 = vfmul_vf_f32m1(vc_10, ALPHA, vl);
      vc_20 = vfmul_vf_f32m1(vc_20, ALPHA, vl);
      vc_30 = vfmul_vf_f32m1(vc_30, ALPHA, vl);
      vse32_v_f32m1(&C[i * N + j], vc_00, tail);
      vse32_v_f32m1(&C[(i + 1) * N + j], vc_10, tail);
      vse32_v_f32m1(&C[(i + 2) * N + j], vc_20, tail);
      vse32_v_f32m1(&C[(i + 3) * N + j], vc_30, tail);
    }
  }
  for (i; i + 1 < M; i += 2) {
    j = 0;
    for (; j + 3 * vl <= N; j += 3 * vl) {
      vfloat32m1_t vc_00 = vle32_v_f32m1(&C[i * N + j], vl);
      vfloat32m1_t vc_01 = vle32_v_f32m1(&C[i * N + j + vl], vl);
      vfloat32m1_t vc_02 = vle32_v_f32m1(&C[i * N + j + 2 * vl], vl);
      vfloat32m1_t vc_10 = vle32_v_f32m1(&C[(i + 1) * N + j], vl);
      vfloat32m1_t vc_11 = vle32_v_f32m1(&C[(i + 1) * N + j + vl], vl);
      vfloat32m1_t vc_12 = vle32_v_f32m1(&C[(i + 1) * N + j + 2 * vl], vl);

      for (k = 0; k < K; ++k) {
        register float a_part_0 = A[i * K + k];
        register float a_part_1 = A[(i + 1) * K + k];
        vfloat32m1_t vb_00 = vle32_v_f32m1(&New_B[j * K + (3 * vl) * k], vl);
        vfloat32m1_t vb_01 =
            vle32_v_f32m1(&New_B[j * K + (3 * vl) * k + vl], vl);
        vfloat32m1_t vb_02 =
            vle32_v_f32m1(&New_B[j * K + (3 * vl) * k + 2 * vl], vl);
        vc_00 = vfmadd_vf_f32m1(vb_00, a_part_0, vc_00, vl);
        vc_01 = vfmadd_vf_f32m1(vb_01, a_part_0, vc_01, vl);
        vc_02 = vfmadd_vf_f32m1(vb_02, a_part_0, vc_02, vl);
        vc_10 = vfmadd_vf_f32m1(vb_00, a_part_1, vc_10, vl);
        vc_11 = vfmadd_vf_f32m1(vb_01, a_part_1, vc_11, vl);
        vc_12 = vfmadd_vf_f32m1(vb_02, a_part_1, vc_12, vl);
      }

      vc_00 = vfmul_vf_f32m1(vc_00, ALPHA, vl);
      vc_01 = vfmul_vf_f32m1(vc_01, ALPHA, vl);
      vc_02 = vfmul_vf_f32m1(vc_02, ALPHA, vl);
      vc_10 = vfmul_vf_f32m1(vc_10, ALPHA, vl);
      vc_11 = vfmul_vf_f32m1(vc_11, ALPHA, vl);
      vc_12 = vfmul_vf_f32m1(vc_12, ALPHA, vl);
      vse32_v_f32m1(&C[i * N + j], vc_00, vl);
      vse32_v_f32m1(&C[i * N + j + vl], vc_01, vl);
      vse32_v_f32m1(&C[i * N + j + 2 * vl], vc_02, vl);
      vse32_v_f32m1(&C[(i + 1) * N + j], vc_10, vl);
      vse32_v_f32m1(&C[(i + 1) * N + j + vl], vc_11, vl);
      vse32_v_f32m1(&C[(i + 1) * N + j + 2 * vl], vc_12, vl);
    }
    for (; j + 2 * vl <= N; j += 2 * vl) {
      vfloat32m1_t vc_00 = vle32_v_f32m1(&C[i * N + j], vl);
      vfloat32m1_t vc_01 = vle32_v_f32m1(&C[i * N + j + vl], vl);
      vfloat32m1_t vc_10 = vle32_v_f32m1(&C[(i + 1) * N + j], vl);
      vfloat32m1_t vc_11 = vle32_v_f32m1(&C[(i + 1) * N + j + vl], vl);

      for (k = 0; k < K; ++k) {
        register float a_part_0 = A[i * K + k];
        register float a_part_1 = A[(i + 1) * K + k];
        vfloat32m1_t vb_00 = vle32_v_f32m1(&New_B[j * K + (2 * vl) * k], vl);
        vfloat32m1_t vb_01 =
            vle32_v_f32m1(&New_B[j * K + (2 * vl) * k + vl], vl);
        vc_00 = vfmadd_vf_f32m1(vb_00, a_part_0, vc_00, vl);
        vc_01 = vfmadd_vf_f32m1(vb_01, a_part_0, vc_01, vl);
        vc_10 = vfmadd_vf_f32m1(vb_00, a_part_1, vc_10, vl);
        vc_11 = vfmadd_vf_f32m1(vb_01, a_part_1, vc_11, vl);
      }
      vc_00 = vfmul_vf_f32m1(vc_00, ALPHA, vl);
      vc_01 = vfmul_vf_f32m1(vc_01, ALPHA, vl);
      vc_10 = vfmul_vf_f32m1(vc_10, ALPHA, vl);
      vc_11 = vfmul_vf_f32m1(vc_11, ALPHA, vl);
      vse32_v_f32m1(&C[i * N + j], vc_00, vl);
      vse32_v_f32m1(&C[i * N + j + vl], vc_01, vl);
      vse32_v_f32m1(&C[(i + 1) * N + j], vc_10, vl);
      vse32_v_f32m1(&C[(i + 1) * N + j + vl], vc_11, vl);
    }
    for (; j + vl <= N; j += vl) {
      vfloat32m1_t vc_00 = vle32_v_f32m1(&C[i * N + j], vl);
      vfloat32m1_t vc_10 = vle32_v_f32m1(&C[(i + 1) * N + j], vl);

      for (k = 0; k < K; ++k) {
        register float a_part_0 = A[i * K + k];
        register float a_part_1 = A[(i + 1) * K + k];
        vfloat32m1_t vb_00 = vle32_v_f32m1(&New_B[j * K + vl * k], vl);
        vc_00 = vfmadd_vf_f32m1(vb_00, a_part_0, vc_00, vl);
        vc_10 = vfmadd_vf_f32m1(vb_00, a_part_1, vc_10, vl);
      }
      vc_00 = vfmul_vf_f32m1(vc_00, ALPHA, vl);
      vc_10 = vfmul_vf_f32m1(vc_10, ALPHA, vl);
      vse32_v_f32m1(&C[i * N + j], vc_00, vl);
      vse32_v_f32m1(&C[(i + 1) * N + j], vc_10, vl);
    }
    if ((tail = N - j)) {
      vfloat32m1_t vc_00 = vle32_v_f32m1(&C[i * N + j], tail);
      vfloat32m1_t vc_10 = vle32_v_f32m1(&C[(i + 1) * N + j], tail);

      for (k = 0; k < K; ++k) {
        register float a_part_0 = A[i * K + k];
        register float a_part_1 = A[(i + 1) * K + k];
        vfloat32m1_t vb_00 = vle32_v_f32m1(&New_B[j * K + tail * k], tail);
        vc_00 = vfmadd_vf_f32m1(vb_00, a_part_0, vc_00, tail);
        vc_10 = vfmadd_vf_f32m1(vb_00, a_part_1, vc_10, tail);
      }
      vc_00 = vfmul_vf_f32m1(vc_00, ALPHA, vl);
      vc_10 = vfmul_vf_f32m1(vc_10, ALPHA, vl);
      vse32_v_f32m1(&C[i * N + j], vc_00, tail);
      vse32_v_f32m1(&C[(i + 1) * N + j], vc_10, tail);
    }
  }
  for (; i < M; ++i) {
    j = 0;
    for (; j + 3 * vl <= N; j += 3 * vl) {
      vfloat32m1_t vc_00 = vle32_v_f32m1(&C[i * N + j], vl);
      vfloat32m1_t vc_01 = vle32_v_f32m1(&C[i * N + j + vl], vl);
      vfloat32m1_t vc_02 = vle32_v_f32m1(&C[i * N + j + 2 * vl], vl);

      for (k = 0; k < K; ++k) {
        register float a_part = ALPHA * A[i * K + k];
        vfloat32m1_t vb_00 = vle32_v_f32m1(&New_B[j * K + (3 * vl) * k], vl);
        vfloat32m1_t vb_01 =
            vle32_v_f32m1(&New_B[j * K + (3 * vl) * k + vl], vl);
        vfloat32m1_t vb_02 =
            vle32_v_f32m1(&New_B[j * K + (3 * vl) * k + 2 * vl], vl);
        vc_00 = vfmadd_vf_f32m1(vb_00, a_part, vc_00, vl);
        vc_01 = vfmadd_vf_f32m1(vb_01, a_part, vc_01, vl);
        vc_02 = vfmadd_vf_f32m1(vb_02, a_part, vc_02, vl);
      }
      vc_00 = vfmul_vf_f32m1(vc_00, ALPHA, vl);
      vc_01 = vfmul_vf_f32m1(vc_01, ALPHA, vl);
      vc_02 = vfmul_vf_f32m1(vc_02, ALPHA, vl);
      vse32_v_f32m1(&C[i * N + j], vc_00, vl);
      vse32_v_f32m1(&C[i * N + j + vl], vc_01, vl);
      vse32_v_f32m1(&C[i * N + j + 2 * vl], vc_02, vl);
    }
    for (; j + 2 * vl <= N; j += 2 * vl) {
      vfloat32m1_t vc_00 = vle32_v_f32m1(&C[i * N + j], vl);
      vfloat32m1_t vc_01 = vle32_v_f32m1(&C[i * N + j + vl], vl);

      for (k = 0; k < K; ++k) {
        register float a_part = ALPHA * A[i * K + k];
        vfloat32m1_t vb_00 = vle32_v_f32m1(&New_B[j * K + (2 * vl) * k], vl);
        vfloat32m1_t vb_01 =
            vle32_v_f32m1(&New_B[j * K + (2 * vl) * k + vl], vl);
        vc_00 = vfmadd_vf_f32m1(vb_00, a_part, vc_00, vl);
        vc_01 = vfmadd_vf_f32m1(vb_01, a_part, vc_01, vl);
      }
      vc_00 = vfmul_vf_f32m1(vc_00, ALPHA, vl);
      vc_01 = vfmul_vf_f32m1(vc_01, ALPHA, vl);
      vse32_v_f32m1(&C[i * N + j], vc_00, vl);
      vse32_v_f32m1(&C[i * N + j + vl], vc_01, vl);
    }
    for (; j + vl <= N; j += vl) {
      vfloat32m1_t vc = vle32_v_f32m1(&C[i * N + j], vl);

      for (k = 0; k < K; ++k) {
        register float a_part = ALPHA * A[i * K + k];
        vfloat32m1_t vb = vle32_v_f32m1(&New_B[j * K + vl * k], vl);
        vc = vfmadd_vf_f32m1(vb, a_part, vc, vl);
      }
      vc = vfmul_vf_f32m1(vc, ALPHA, vl);
      vse32_v_f32m1(&C[i * N + j], vc, vl);
    }
    if ((tail = N - j)) {
      vfloat32m1_t vc = vle32_v_f32m1(&C[i * N + j], tail);

      for (k = 0; k < K; ++k) {
        register float a_part = ALPHA * A[i * K + k];
        vfloat32m1_t vb = vle32_v_f32m1(&New_B[j * K + tail * k], tail);
        vc = vfmadd_vf_f32m1(vb, a_part, vc, tail);
      }
      vc = vfmul_vf_f32m1(vc, ALPHA, vl);
      vse32_v_f32m1(&C[i * N + j], vc, tail);
    }
  }

  free(New_B);
}

/**
 * @brief 根据l2 cache计算分块矩阵的大小。
 * @param[in]     M     矩阵A的行数
 * @param[in]     N     矩阵B的列数
 * @param[in]     K     矩阵A的列数和矩阵B的行数
 * @param[out]    T_M   获取矩阵A分块行数的指针
 * @param[out]    T_N   获取矩阵B分块列数的指针
 * @param[out]    T_K   获取矩阵A分块列数/矩阵B分块行数的指针
 * @param[in]     vl    向量长度
 * @return void
 */
static inline void count_tile_mnk_1(const int M, const int N, const int K,
                                    int *T_M, int *T_N, int *T_K, int vl) {
  int tile_size = (int)sqrtf((float)total_l2_cache_size / 3);

  int TILE_M = max(4, tile_size / 4 * 4);
  int TILE_N = max(vl, tile_size / vl * vl);
  int TILE_K = max(4, tile_size / 4 * 4);

  int nn_K = (K + TILE_K - 1) / TILE_K;
  TILE_K = min(TILE_K, ((K + nn_K - 1) / nn_K + 3) / 4 * 4);

  if (nn_K == 1) {
    tile_size = (int)((float)total_l2_cache_size / 2 / TILE_K);
    TILE_M = max(8, tile_size / 8 * 8);
    TILE_N = max(vl, tile_size / vl * vl);
  }

  *T_M = TILE_M;
  *T_N = TILE_N;
  *T_K = TILE_K;
}


/**
 * @brief 使用顺序读跨步写的方式将矩阵B打包到New_B中，服务于分块算法。
 * @param[in]     B     矩阵B的指针
 * @param[out]    New_B 输出矩阵New_B的指针
 * @param[in]     K     矩阵B的行数
 * @param[in]     N     矩阵B的列数
 * @param[in]     vl    向量长度
 * @return void
 */
static void pack_B_MNK(const float *B, float *New_B, int K, int N, int T_K,
                       int T_N, int nn_K, int nn_N, int vl) {
  int j, k, tail;
  int jj, kk;

  int tail_K = K - T_K * (nn_K - 1);
  int tail_N = N - T_N * (nn_N - 1);

  for (k = 0; k < K; ++k) {
    for (j = 0; j + vl <= N; j += vl) {
      int tile_index_K =
          (j / T_N) * ((nn_K - 1) * (T_K * T_N) + (tail_K * T_N));
      int tile_index_N = N - (j / T_N * T_N) < T_N ? (k / T_K) * tail_N * T_K
                                                   : (k / T_K) * T_K * T_N;
      int tile_index = tile_index_K + tile_index_N;
      int new_b_index = tile_index +
                        (j % T_N) * (min(T_K, K - (k / T_K * T_K))) +
                        vl * (k % T_K);
      int b_index = k * N + j;
      vfloat32m1_t vb = vle32_v_f32m1(&B[b_index], vl);
      vse32_v_f32m1(&New_B[new_b_index], vb, vl);
    }
    if ((tail = N - j)) {
      int tile_index_K =
          (j / T_N) * ((nn_K - 1) * (T_K * T_N) + (tail_K * T_N));
      int tile_index_N = N - (j / T_N * T_N) < T_N ? (k / T_K) * tail_N * T_K
                                                   : (k / T_K) * T_K * T_N;
      int tile_index = tile_index_K + tile_index_N;
      int new_b_index = tile_index +
                        (j % T_N) * (min(T_K, K - (k / T_K * T_K))) +
                        tail * (k % T_K);
      int b_index = k * N + j;
      vfloat32m1_t vb = vle32_v_f32m1(&B[b_index], tail);
      vse32_v_f32m1(&New_B[new_b_index], vb, tail);
    }
  }
}

/**
 * @brief 分块版本的GEMM实现。
 * @param[in]     M     矩阵A的行数
 * @param[in]     N     矩阵B的列数
 * @param[in]     K     矩阵A的列数和矩阵B的行数
 * @param[in]     ALPHA 标量系数
 * @param[in]     A     矩阵A的指针
 * @param[in]     B     矩阵B的指针
 * @param[in,out] C     输出矩阵C的指针
 * @return void
 */
void gemm_nn_rvv_MNK(int M, int N, int K, float ALPHA, const float *A,
                     const float *B, float *C) {
  // 计算每个向量操作可以处理的元素数
  size_t vl = get_max_vlen_e32m1();
  int tail = 0;
  int i, j, k;
  int ii, jj, kk;

  int T_M, T_N, T_K;
  count_tile_mnk_1(M, N, K, &T_M, &T_N, &T_K, vl);


  int nn_K = (K + T_K - 1) / T_K;
  int nn_M = (M + T_M - 1) / T_M;
  int nn_N = (N + T_N - 1) / T_N;

  int tail_K = K - T_K * (nn_K - 1);
  int tail_N = N - T_N * (nn_N - 1);

  float *New_B = (float *)aligned_alloc(64, K * N * sizeof(float));
  pack_B_MNK(B, New_B, K, N, T_K, T_N, nn_K, nn_N, vl);

  vfloat32m1_t vc;
  vfloat32m1_t vb;
  for (ii = 0; ii < M; ii += T_M) {
    int max_ii = min(T_M, M - ii);
    for (jj = 0; jj < N; jj += T_N) {
      int max_jj = min(T_N, N - jj);
      for (kk = 0; kk < K; kk += T_K) {
        int max_kk = min(T_K, K - kk);
        int new_b_index_1 = ((jj / T_N) * ((nn_K - 1) * (T_K * T_N) + (tail_K * T_N))) + (N - (jj / T_N * T_N) < T_N
                                     ? (kk / T_K) * tail_N * T_K
                                     : (kk / T_K) * T_K * T_N);
        for (i = ii; i < ii + max_ii; ++i) {
          int new_b_index = new_b_index_1;
          for (j = jj; j + vl <= jj + max_jj; j += vl) {
            vc = vle32_v_f32m1(&C[i * N + j], vl);
            for (k = kk; k < kk + max_kk; ++k) {
              float a_part = ALPHA * A[i * K + k];
              vb = vle32_v_f32m1(&New_B[new_b_index], vl);
              vc = vfmadd_vf_f32m1(vb, a_part, vc, vl);
              new_b_index += vl;
            }
            vse32_v_f32m1(&C[i * N + j], vc, vl);
          }
          if ((tail = jj + max_jj - j)) {
            vc = vle32_v_f32m1(&C[i * N + j], tail);
            for (k = kk; k < kk + max_kk; ++k) {
              float a_part = ALPHA * A[i * K + k];
              vb = vle32_v_f32m1(&New_B[new_b_index], tail);
              vc = vfmadd_vf_f32m1(vb, a_part, vc, tail);
              new_b_index += tail;
            }
            vse32_v_f32m1(&C[i * N + j], vc, tail);
          }
        }
      }
    }
  }

  free(New_B);
}
